####### Framing Expriment of securitization Turkey 2024 coded by Shingo Hamanaka 20240806 #########
rm(list=ls())

library(tidyverse)
library(haven)
library(stargazer)
library(BalanceR)
library(ggpubr)

data1 <- read.csv("~/Turkey_securitization2024_num.csv")


#### outcome：数字が大きいほど支持しないので逆順に ####
data2<-
  data1 %>%
  mutate(Q44R = - Q44 + 5)
data3<-
  data2 %>%
  mutate(Q45R = - Q45 + 5)
data4<-
  data3 %>%
  mutate(Q47R = - Q47 + 5)
data5<-
  data4 %>%
  mutate(Q48R = - Q48 + 5)

data5$Q44R[is.na(data5$Q44R)]<-0
data5$Q45R[is.na(data5$Q45R)]<-0
data5$Q47R[is.na(data5$Q47R)]<-0
data5$Q48R[is.na(data5$Q48R)]<-0

data5$outcome<-data5$Q44R+data5$Q45R+data5$Q47R+data5$Q48R
data5$outcome[data5$outcome==0]<-NA

#### treatment ###
data5<-data5 %>% 
  mutate(
    treatment =
      case_when(Q44R>=1 ~ "Picture",
                Q45R>=1 ~ "Securitization",
                Q47R>=1 ~ "Picture+Securitization",
                Q48R>=1 ~ "Control"
      )
  )
data5$treatment <- as.factor(data5$treatment)  ###因子化：1は写真、2は安全保障化言説、3は写真と言説、4は統制群

### Balance Check ####
library(BalanceR)
balance_chk <- BalanceR(data = data5, group = treatment,
                        cov = c(Q2, Q3, Q4, Q6))
print(balance_chk, digits = 3)
summary(balance_chk, digits = 5)
plot(balance_chk, point.size = 5, text.size = 18)



#### 分散分析 ###
pairwise.t.test(data5$outcome, data5$treatment, p.adj="none")

TukeyHSD(aov(data5$outcome~data5$treatment))

library(multcomp)


res<- lm(outcome~treatment,data=data5)
tukey_res<- glht(res, linfct=mcp(treatment="Tukey"))
summary(tukey_res)

pairwise.t.test(data5$outcome, data5$treatment, p.adj="bonf")

pairwise.t.test(data5$outcome, data5$treatment, p.adj="BH")

kruskal.test(outcome~treatment,data=data5)

#### graph ####

data5$Q44R[data5$Q44R==0]<-NA
data5$Q45R[data5$Q45R==0]<-NA
data5$Q47R[data5$Q47R==0]<-NA
data5$Q48R[data5$Q48R==0]<-NA

summary(data2$Q44R)
summary(data3$Q45R)
summary(data4$Q47R)
summary(data5$Q48R)

framing_result1 <- data5 %>% 
  summarise(pic_mean = mean(Q44R, na.rm = TRUE), pic_sd = sd(Q44R, na.rm = TRUE),
            word_mean = mean(Q45R, na.rm = TRUE), word_sd = sd(Q45R, na.rm = TRUE),
            pic_word_mean = mean(Q47R, na.rm = TRUE), pic_word_sd = sd(Q47R, na.rm=TRUE),
            cntl_mean = mean(Q48R, na.rm=TRUE), cntl_sd = sd(Q48R, na.rm=TRUE))

sd(data5$Q44R, na.rm = TRUE) / sqrt(length(data5$Q44R))
sd(data5$Q45R, na.rm = TRUE) / sqrt(length(data5$Q45R))
sd(data5$Q47R, na.rm = TRUE) / sqrt(length(data5$Q47R))
sd(data5$Q48R, na.rm = TRUE) / sqrt(length(data5$Q48R))


p1 <- data.frame(
  Experiment = c("picture", "securitization", "pic+securi", "control"),
  Mean = c(2.712, 2.816,  2.803, 2.704),
  SE = c(0.02863036, 0.02703427, 0.02778221, 0.02709348)
)
ggplot(p1, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean + SE, ymin = Mean - SE), 
                width = 0.2, position = position_dodge(.9)) +
  labs(title = "Framing experiment of securitization") +
  theme_bw()


#### treatment 2 ###
data5<-data5 %>% 
  mutate(
    treatment2 =
      case_when(Q44R>=1 ~ "Control",
                Q45R>=1 ~ "Securitization",
                Q47R>=1 ~ "Securitization",
                Q48R>=1 ~ "Control"
      )
  )
data5$treatment2 <- as.factor(data5$treatment2)  ###因子化

### Balance Check ####
library(BalanceR)
balance_chk <- BalanceR(data = data5, group = treatment2,
                        cov = c(Q2, Q3, Q4, Q6))
print(balance_chk, digits = 3)
summary(balance_chk, digits = 5)
plot(balance_chk, point.size = 5, text.size = 18)


#### 分散分析 ###
pairwise.t.test(data5$outcome, data5$treatment2, p.adj="none")

TukeyHSD(aov(data5$outcome~data5$treatment2))

res2<- lm(outcome~treatment2,data=data5)
tukey_res2<- glht(res2, linfct=mcp(treatment2="Tukey"))
summary(tukey_res2)

pairwise.t.test(data5$outcome, data5$treatment2, p.adj="bonf")

pairwise.t.test(data5$outcome, data5$treatment2, p.adj="BH")


wilcox.test(outcome~treatment2,data=data5)

